Genomic identification and characterization of prophages associated with Citrobacter freundii strains

Citrobacter freundii belongs to the Enterobacteriaceae family of facultative anaerobic Gram-negative bacteria [1]. They are soil-dwelling microbes in water, sewage, food, animal, and human digestive systems, sometimes responsible for community-acquired and nosocomial infections [2, 3]. C. freundii can cause a life-threatening infection that can be turned into sepsis and has also been associated with brain abscesses and newborn meningitis [4, 5]. The virulence factors of C. freundii include adhesins, endotoxin, hemolysin, and siderophores, which help it to colonize and invade host tissues [6]. Furthermore, the resistance of C. freundii to various antibiotics is also a growing concern in healthcare settings [5, 7].
Bacteriophages are viruses that infect bacteria and are abundant in the biosphere [8]. Lytic or virulent and lysogenic or temperate are the two types of bacteriophages [9]. The tailed phages are the most prominent group of viruses that affect bacteria and are the most abundant in the human gut belonging to the order Caudovirales [10]. Although the initial roles of phage genes are to infect, survive, and reproduce in the host bacteria successfully, several pathogenic and antibiotic-resistant genes are encoded by the phage genomes. Therefore, understanding the biology of the phages has several implementations in selecting phages as antibacterial agents, the evolution of bacterial strains, and typing for epidemiological studies [10]. Previously, characterization and comparative analysis of prophages present in several pathogenic strains of Burkholderia cenocepacia [11], Helicobacter pylori [12], Streptococcus mutans [13], Acinetobacter baumannii [14, 15], Acinetobacter pittii [16], Klebsiella pneumoniae [17, 18], Desulfovibrio [19], Staphylococcus aureus [20] have been performed to understand the roles of prophages in the evolution, fitness, and pathogenicity of the host bacteria. 
This study elucidated the genomic characterization of 574 prophage genomes from C. freundii. Moreover, various comparative bioinformatics tools were used to study the phylogenomic, evolutionary, and virulence profile analysis of the intact prophages.

Identification and classification of prophages sequences

The complete chromosomal sequences of 144 C. freundii strains, submitted to the National Center of Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/genome/browse/#!/prokaryotes/2850/), were retrieved until November 2022 as described in [16, 17]. Briefly, from the database C. freundii genome sequences were searched. Later, only complete C. freundii genome sequences were downloaded.  Plasmid genomes were not included in the analysis. The genome size of each strain was obtained from GenBank. All the strains were screened for prophage-like elements utilizing the web server Phage Search Tool Enhanced Release (PHASTER) [21]. PHASTER was used to identify prophages as intact, questionable, or incomplete based on the length of the phage and the presence of phage-derived genes. The integrity of prophage-associated genes like integrase, endolysin, and terminase genes and the 'att' sites (attL/attR) in intact phages were manually checked for all intact prophages to reduce the likelihood of false positives. Prophage families of intact prophage were characterized by the Virfam server based on their structural compositions [22]. The recognition of head-neck-tail modules was used in Virfam for the classification of the viruses.

Prophage-associated virulence factors and antibiotic-resistance genes

Annotation of the protein sequences present in the intact prophages was performed by Rapid Annotation using System Technology (RAST) version 2.0 [23]. In Microsoft Excel, a heatmap was generated by using conditional formatting of the cells based on toxin proteins present in the intact prophages. For the prediction of antibiotic resistance genes, genomic sequences of prophages were screened by using the resistance gene identifier (RGI) tool from the comprehensive antibiotic resistance database (CARD) [24] and the ABRicate tool from Galaxy Australia [25]. The default parameters were used in all analyses conducted.

Phylogenomic analysis of intact prophage sequences

The phylogenomic study of all intact prophages was performed in Molecular Evolutionary Genetics Analysis (MEGA X) [26] using Neighbouring Joining methods with 1000 bootstrap replications. The tree was visualized by Interactive Tree of Life (iTOL) [27].

Spatial distribution of intact prophages

The spatial distribution of intact prophages inside bacterial chromosomes was performed following the protocol described in [16] using the R package ggplot2 [28]. To generate a density plot of the intact prophages, lengthwise bacterial chromosome sequences of the geom_density () function of ggplot2 were used.

Prevalence of prophages in the C. freundii genomes

PHASTER predicted 1178 prophages from 144 C. freundii strains isolated from various sources of different countries. One strain lacks intact prophage (C. freundii CitB, Accession No: CP048416). PHASTER predicted 604 incomplete and questionable prophages and 574 intact prophages (average 3.98±1.96) (Figure 1A, Supplementary Table 1). Each chromosome, on average, had 5.2% of the total prophage genome present in C. freundii. For further analysis, we considered only intact prophages. The intact prophage sequence length ranged from 6.74 kb to 115.15 kb (Figure 1B). However, around 52% of prophages had genome sizes ranging from 41-60 kb, only one prophage ranged from 81-100 kb, and 4 prophages had genome sizes between 101-120 kb (Figure 1B). Most frequently, 3 intact prophages were identified in 38 (26.39% of 144) C. freundii. A maximum of 9 intact prophages were detected in three C. freundii (2.1%) (Figure 1C).

Classification of intact prophages in C. freundii strains

The classification of 574 intact prophages in various viral families was performed according to the proteins of head and tail connectivity, a parameter used by the VirFam tool. The majority of intact prophages (48.16%) were predicted to belong to the Myoviridae family, followed by the Siphoviridae (42%), Podoviridae (4.67%), and Inoviridae (0.17%) families. The server predicted 31(5%) unclassified prophages. The average prophage genome sizes of the Podoviridae and Siphoviridae were slightly higher than Myoviridae (Figure 2). In contrast, smaller genome sizes and lower GC contents were observed in the prophage in Inoviridae than others. This observation was due to only one prophage identified from the Inoviridae family.

Later, a density plot for 574 intact prophages present in C. freundii chromosomes was generated. This density plot demonstrated a significant prophage density of around 54900000 bp and 5100000 bp, as shown in Figure 3.

Presence of virulence and toxin proteins within prophages

Protein-coding sequences of the intact prophages were analyzed for the presence of virulence-associated proteins. Out of 574 intact prophages, 152 (26.43%) harbored virulence genes. However, out of the 152 intact prophages, 118 contained MsgA (Macrophage-specific gunine nucleotide Exchange Factor-like) protein (Supplementary Table 2), 10 prophages were predicted to have YkfI proteins, 11  had RTX (Repeats-in-toxin) exo-proteins, 7 had HigB (one of several type III secreted effector proteins produced by E. coli, and in P. aeruginosa, HigB modulates the expression of T3SS genes and biofilm formation) proteins, 4 had YeeV proteins, and 1 prophage had Hcp (Haemolysin Co-regulated Protein) and another had HipA (high persistence A)  proteins (Figure 4).

Phylogenomic analysis of the intact prophages

The whole genome sequences of 574 prophages were used to analyze the phylogenomic relationship. The results revealed that these intact prophages were classified into three major clades (Figure 5). Clade 1, clade 2, and clade 3 have 105, 155, and 314 prophage genomes, respectively. Besides, 22.86% of prophages in clade 1, 21.29% in clade 2, and 30.25% in clade 3 have virulent genes (Supplementary Table 3).

Numerous bacterial isolates have prophages present, and these prophages may have contributed to the diversity of bacteria. Prophages can contribute to a bacterial strain's genetic distinctiveness by carrying different virulent components. Among other organisms, this phenomenon has been noted in Lactobacillus [29], Pseudomonas [30], Salmonella genus [10, 31], Streptococcus pyogenes [32], Streptococcus mutans [13],  K. pneumoniae [17], Paenibacillus larvae [33]. Previous studies have shown that strain-specific differences are caused mainly by prophage-associated chromosomal rearrangements and deletions [17, 30, 31, 34]. Therefore, the emergence of pathogenic variants may be influenced by temperate phages carrying foreign DNA. The production of toxins and disease-modifying substances encoded by phages and phage lysogenic conversion genes can change the host phenotype and affect the virulence of the bacterial strain [35].

In this study, 1178 prophages were predicted in 144 C. freundii strains, with an average of 5.2% of the host chromosome content prophage genome, which is lesser than other bacterial genomes (10-20%) [35, 36]. On the other hand, a more significant fraction of prophages in C. freundii were found to be defective (51.27%). The incomplete and questionable prophages frequently lack essential phage functions. Such deficiency might be caused by the selective pressure of the host to inactivate mutations on the newly acquired phages. Consequently, the causative genes of infection were deleted or mutated, providing multiple benefits to the host bacterium to adjust to hostile environmental conditions [37]. Therefore, our further analysis focused on intact prophages characterized by an integrase domain, transcriptional regulator, or at least one of the structural domains attributed to either capsid, tailed, or portal proteins. The prophage GC content is 50.75%, lower than the average GC content of 51.85% in the C. freundii chromosomes. The intact prophage sizes in C. freundii range from 6.1 kb to 399.9 kb. In contrast, another member of the Enterobacteriaceae family had the prophage sizes in K. pneumoniae range from 29.7 kb to 50.6 kb [18], 6.4 kb to 358 kb for Salmonella [31], and 8.1 kb to  31.9 kb for Shigella [38]. There seems to be a massive difference in the prophage sizes of the members of the Enterobacteriaceae family. Similar prophages in different host isolates indicate their horizontal movement events and significance in genomic plasticity or evolution (Supplementary Table 4) [39]. C. freundii strains were shown to have multiple phages ranging from 10 to 17, belonging to different families and clades in different groups. Multiple prophage-like elements suggest that poly-lysogeny may be shared in many bacterial pathogens, including C. freundii.

The viral families are the basic unit of the biological classification of viruses, and these are characterized mainly based on physiochemical properties, genome structures, morphology, and gene expression [22]. Among the 574 intact prophages studied, the In silico approaches revealed Myoviridae as the most represented family (48.16%), followed by Siphoviridae (42%), Podoviridae (4.67%) and Inoviridae (0.17%). Similar to our studies, another investigation focusing on the prophages of K. pneumoniae revealed that 62 prophages in Myoviridae (59.6%), followed by  40 prophages in  Siphoviridae (38.5%) and, 2 prophages in Podoviridae (1.9%) could be classified [18]. Furthermore,  the prophages in A. baumannii could be divided into Myoviridae (25.69%), Siphoviridae (31.20%),  Podoviridae (10.09%), and unclassified (33.02%) [40]. In C. freundii, Myoviridae is typically the most prominent phage family and has below-average genome sizes (42.25 kb). In comparison, Siphoviridae (45.62 kb) and Podoviridae (46.25 kb) had sizes slightly above the average (42.97 kb) (Figure 2). A phylogenomic tree was generated and grouped into three major clades to understand the identified diversity of intact prophages (Figure 5, Supplementary Table 3). All the clades were revealed to be a mixed sub-cluster with higher diversity, comprising prophages from Myoviridae, Siphoviridae, Podoviridae, Inoviridae, and unclassified families based on the program used in the study.

Among the 574 intact prophages, the most dominating are Klebsiella phage phiKO2 (9.93%) and Salmonella phage 118970_sal3 (9.58%) (Supplementary Table 4). These two phages are present in several C. freundii strains isolated from several locations. The 574 intact prophages studied showed that some shared more than 95% of their genome identities, indicating solid evolutionary relationships. We discovered 65 phages that shared more than 70% similarities with the 574 intact prophages regarding identity when compared to the studied prophage genomes against NCBI databases (Supplementary Table 4).

Multiple studies have shown the relationship between the pathogen's virulence and the acquisition of prophages leading to numerous infectious diseases. Some prophages carry genes that can change the characteristics of the host, ranging from increased host fitness to virulence [41-43]. Even defective prophages are thought to be capable of transporting virulence factors [44]. In our 152-virulence gene containing intact prophages analysis, msgA was predominantly distributed. (Supplementary Table 2). The MsgA exhibits a distinctive expression pattern during macrophage infection that contributes to yeast cell morphogenesis variation [45]. Additionally,  MsgA has been suggested to play a role in the biofilm formation and antibiotic resistance in Bacillus subtilis and S. aureus [46, 47]. Therefore, this protein may contribute to biofilm formation and antibiotic resistance in C. freundii.

Some C. freundii prophages have been found to have T1SS-secreted agglutinin RTX protein with diverse biological functions like virulence, biofilm formation, and resistance to antibiotics [48]. But in Citrobacter spp. the activities of T1SS-secreted agglutinin RTX proteins are still unknown, although the prevalence of RTX protein was the second highest (~6.57%) after the MsgA virulence factor identified in the current study. A previous study showed that RTX-containing proteins often have hemolytic or hemagglutinating activity in B. cenocepacia J2315 and identified four hemolysin-type Ca²⁺-binding domains at the carboxy-terminal of the protein sequence that is characteristic of RTX-related proteins [49]. Therefore, some prophages confer new virulence factors and/or antibiotic resistance genes used for bacterial pathogenesis and beneficial traits like increased fitness to the host [44, 50].

Furthermore, we have found a gene for the Hcp family protein in one of the intact prophages (Figure 4). The Hcp family of proteins is secreted through the Type 6 Secretion System (T6SS), which plays a crucial role in the pathogenesis of bacterial infections and contributes to the virulence of the invading bacteria, including E. coli, Vibrio, Aeromonas, and Pseudomonas [51-53]. According to previous research, the T6SS genomic island in  C. freundii CF74 strain modulates the flagellar system, improves motility, thus is involved in adhesion to host cells, and produces cytotoxicity in host cells [54]. This study detected the toxin gene hipA in one of the intact prophages from C. freundii strain 705SK3. HipA toxin is a part of the HipA-HipB Type II Toxin Antitoxin (TA) system, where HipA is involved in the inhibition of Glutamyl-tRNA synthetase, thereby ceasing protein synthesis and leading to the death of the bacteria [53, 55]. HipA-HipB Type II, TA system, was reported in the IncA/C1 Plasmid pCf587 of C. freundii 33587 [56].

The heat map analysis revealed the presence of another three different proteins, Ykfl, HigB, and YeeV, among the remaining intact prophages associated with virulence factors. One study conducted using E. coli YeeV toxin reported that YeeV interacts with two major cytoskeleton proteins, FtsZ (a tubulin-like protein) and MerB (a prokaryotic actin-like protein), to inhibit the activities of both GTPase and GTP-dependent polymerization of FtsZ as well as the ATP-dependent polymerization of MerB. Consequently, it inhibits cell division and changes the morphology and lysis of E. coli [57]. Another study reported that prolonged overproduction of YkfI causes complete cell death within seconds or minutes by interacting with the cell division protein FtsZ in cryptic prophages of E. coli K-12 [58]. It is also evident from the same study that toxin YkfI constitutes a toxin-antitoxin (TA) pair with type IV antitoxin YafW. Such antitoxin can only partially block the toxicity of YkfI. However, the interaction of  YkfI with FtsZ did not cause cell death in Pseudomonas strains [58]. Several studies reported the role of toxin HigB, a member of the ribosome-dependent RelE toxin family, in association with Type II antitoxin HigA in regulating pathogenicity in many opportunistic pathogenic bacteria, including  E.coli K-12, P. aeruginosa PA14, Vibrio cholerae,  Streptococcus pneumoniae [55, 59, 60]. Therefore, further research is needed to determine whether C. freundii prophages can benefit the host's evolutionary fitness due to these virulence factors.

To summarize, this study successfully identified the intact prophages in all completely sequenced strains of C. freundii that were previously unknown. This study also investigated the genomes of these prophages to determine their potential impact on the host strains. C. freundii, an opportunistic pathogen, harbors prophages that could affect its virulence and cause genetic rearrangements. Nearly half of the identified prophages were incomplete, suggesting that C. freundii is actively integrating and modifying its prophages. The intact prophages may actively influence the bacterial genome and contribute to the pathogenicity of C. freundii by encoding significant virulence genes. Many prophages may be crucial in bacterial evolution, adaptability, and physiology. Furthermore, conducting prophage induction experiments would provide insights into host-phage interactions. The research findings also indicate the presence of multiple prophage-like elements, known as poly-lysogeny, which can convert non-pathogenic bacteria into pathogenic ones capable of producing harmful toxins. Additionally, the results demonstrate that C. freundii and other members of the Enterobacteriaceae family share similar patterns in terms of prophage distribution and the prevalence of common prophages.

Authors would like to acknowledge the Department of Biochemistry and Microbiology, North South University for their technical support for the project.

Conceptualization: SRS; Data curation: SIM, ABL, AA, SHO; Formal analysis: SIM, AR, ABL, IJ; Methodology: SRS, ABL, SIM, IJ; Writing – original draft: IJ, SIM, SI, SRS; and Writing – review & editing: SRS, IJ, SI. All authors have read and approved the final version of the manuscript.

There is no conflict of interest among the authors.

Supplementary Table 1. Characterization of the prophage genomes in C. freundii genomes, Supplementary Table 2. List of Toxin and virulence proteins, Supplementary Table 3. Distribution of the prophages in different clades, and Supplementary Table 4. Representative prophages present in the genomes of C. freundii (Supplementary Materials).